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Abstract. Analyzing in detail the first corrections to the scaling hypothesis, we develop 
accelerated methods for the determination of critical points from finite size data. The output 
of these procedures is sequences of pseudo-critical points which rapidly converge towards the 
true critical points. In fact, more rapidly than previously existing methods such as the 
phenomenological renormalization group (PRG) approach. Our methods are valid in any spatial 
dimensionality and both for quantum or classical statistical systems. Having fast converging 
sequences at our disposal, allows us to draw conclusions on the basis of shorter system sizes. 
This can be extremely important in particularly hard cases such as two-dimensional quantum 
systems with frustrations or in Monte Carlo simulations when the sign problem occurs. We test 
the effectiveness of our methods both analytically on the basis of the one-dimensional XY model, 
and numerically at phase transitions occurring in non integrable spin models. In particular, we 
show how a new Homogeneity Condition Method is able to locate the onset of the Berezinskii- 
Kosterlitz-Thouless transition simply making use of ground-state quantities on relatively small 
systems. 
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1. Introduction 



In statistical mechanics the exact evaluation of equilibrium properties, except for non-interacting 
systems, can be worked out only in very few notable integrable models [6||. Thanks to the increasing 
power of computer resources, numerical simulations on finite size systems are becoming a more and 
more popular. Among them we quote exact diagonaHzations (ED) algorithms, density matrix 
renormalization group (DMRG), and (quantum) Monte Carlo methods. The first question one asks 
in this context is the determination of critical points in order to classify the phase diagram into 
different phases. Once critical points are found, one can ask what is the nature of the critical point 
by computing critical exponents. Though statistical mechanics was born more than a hundred years 
ago, and its theoretical understanding is now on a sound basis after the scaling hypothesis and the 
theory of renormalization group, until now only somehow phenomenological methods are available 
for the determination of critical points from finite size data. In this article we carefully analyze the 
leading corrections of equilibrium averages beyond the usual scaling hypothesis. In this way we are 
able to find accelerated methods for the determination of critical points. The major advantage of 
having faster converging methods is that one can content oneself with simulations on shorter sizes. 
Such simulations are those which can be obtained with higher accuracy in favorable cases, but are 
the only one available in the most demanding situations given for instance by frustrated quantum 
systems in two dimensions. 

Traditionally classical (some would say Euclidean) statistical systems, have been studied by 
signaling out an order parameter (when present). An example can be the magnetization in a 
magnetic system (M), which, in the thermodynamic limit, is non zero only in the ordered phase, 
and so allows to locate precisely the critical point. This is only true in the thermodynamic limit, 
instead at finite sizes the situation is less clear, and typically the magnetization will be different 
from zero everywhere, presenting "finite size tails" in the paramagnetic region. A possible, widely 
used, strategy is that of concentrating on susceptibilities, i.e. derivative of the magnetization. It 
can happen that the susceptibility diverges at the critical point in the infinite volume limit. At 
finite size singularity are smoothed and the susceptibility only displays a maximum at a certain 
point which is shifted respect to the actual critical point. Of course when one considers larger and 
larger sizes, the pseudo-critical point given by the maximum of the susceptibility, shifts towards 
the true critical point. The convergence rate is algebraic in the system size L and is of the form 
defining the so called shift exponent 9. This method of course is limited to those cases where 
a susceptibility diverge. A somewhat more elaborated method is that of the Binder's cumulant ■ 

In this approach one considers the fourth order reduced cumulant J/^ = 1 — ((AM)**)/ 3((AM)^)^ 

where AM = M — (M), in a finite system of linear size L. Since by the scaling hypothesis Ul 
should depend only on the ratio L/^ where ^ is the correlation length, at the critical point ^ — > oo 
and U L should be independent of L. The finite size pseudo critical point is then obtained by solving 
SlUl = (here and in the following, at finite size, derivatives with respect to L must be interpreted 
as finite differences) . Graphically this amounts to finding the crossing points of the curves Ul at 
successive values of L. 

In quantum statistical systems the method of choice is given by the so called phenomenological 
renormalization group (PRG) [1, Continuous phase transitions in quantum mechanics are 

characterized by closure of the first gap A/,. Scaling hypothesis then asserts that the gap must be 
a function of ^ hence A^ ~ C"'' which defines the dynamical exponent C,. Since at the critical point 
the correlation length will be of the order of the system size, the quantity L'' A^ will be independent 
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of L. The PRG equation is then given by Ol [L'^Al] = 0. 

All these methods produce a sequence of pseudo critical points, say r£ which scale to the critical 
point according to an algebraic law T£ = Tc + O (L^^) governed by the shift exponent 9. The 
general belief is that the shift exponent 9 is given by the inverse of the correlation length exponent 
ly. Although, to our knowledge, no systematic study of shift exponents has been done, it was 



already pointed out in that the relation 9 — l/v is not always valid, and 9 depends, among 
other factors, on the boundary conditions. In particular it should be clear that 9 depends on the 
specific way in which the pseudo-critical point has been obtained. 

It should be reminded here that in recent times new tools for the characterization of quantum 
critical points, inherited from the theory of quantum information, have arisen Among these 
the so called fidelity approach to phase transition consists of looking for the maximum of the so 
called quantum Fisher information [9, 22]. This method is the best possible one if the only available 
resource we admit is measuring observables. As we will see we can do much more, i.e. compare 
different observables, compute finite difference and so on. In any case it was shown in |23l] that, 
except for accidental cancellations, the shift exponent obtained with this method is 9 = 1/v. 

The lesson we should bear from this short resume, is twofold. On the one hand the present methods 
for the location of critical point are somehow phenomenologically obtained assuming/imposing the 
scaling hypothesis. On the other hand scaling hypothesis alone would imply that finite size pseudo- 
critical points actually coincide with critical points. To learn about the shift exponents and to 
develop faster converging sequences one must go beyond the usual scaling hypothesis. 



2. Finite size scaling beyond scaling hypothesis 

To be concrete we will deal with quantum systems in d spatial dimension at zero temperature 
endowed with periodic boundary conditions (PBC) to get rid of border effects. Our results can 
however be easily generalized to finite temperature (in which case the role of the energy will be 
played by the free energy) and/or classical statistical mechanics in D dimension (and use the 
dimensional crossover rule D ^ d + (). The Hamiltonian can be written as 

n^no+gV = ni + {g-gc)V. 

At the critical point gc the model is given by the Hamiltonian Hi = Ho+gcV which, by hypothesis, 
is critical and so defines a scale invariant theory. The perturbation V can then be expanded in 
terms of the scale-invariant operators of the theory Hi, this is the so called hypothesis of local 
algebra [l3|. These operators are classified into relevant /irrelevant in case their scaHng dimensions 
are smaller/larger than d + ( , the case d + being marginal. Since, for g ^ gc the model is not 
critical, there must be at least one relevant component in V which is responsible for driving the 
system away from criticality. This simple observation motivates the fact that, at least in principle, 
it should be possible to estimate the critical point, only considering the finite size sequence of the 
perturbation operator hh = (V) jL"^. Such a strategy has been advocated in where the so called 
finite size crossing method (FSCM) has been developed. It was found that a sequence of pseudo- 
critical points g*j^ could be given by the simple FSCM equation djjah (ffl,) = 0. As we will see later 
this method results in a shift exponent twice as large as the "standard" one: 9pscM = ^/v. 
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The energy density e = (H) / L"^ can be splitted into a regular and a singular part: e — Crog + Csing- 
In the ofF-critical regime L ^ £^ the singular part behaves as 

esing - (1) 

where t = g — gc- In the most general case, the coefficient in front of Eq. |(T]) can depend on the 
side the transition is approached, i.e. Csmg — A± \t\ ~" according to whether t ^ 0"^. ScaHng 
hypothesis asserts that, in a finite system of linear size L, thermodynamic averages depend only on 
the ratio L/£^, where, close to the transition point, ^ ~ \t\ ^ . Since for finite systems averages must 
be analytic in g, the singular part of the energy must effectively depend on the scaling variable 
z = tL^/^ . So we arrive at 

esi„g = i-(''+'^)<I>oW + ---, (2) 

where the dots represent corrections arising from less singular contributions. Here $0 [z) is a 
universal scaHng function that for z ±cx3 (the off-critical regime) must behave as $0 [z) ~ 
A_^z^'i+Qu order to recover Eq. Jl]). Comparing with Eq. (U) one finds 2 — a— {d + C,)v. Instead, 
for L ^ ^ we are in the critical regime and <i>o(-z) behaves as an analytic function. The only 
requirement is that, at finite L, $0 must be such that e is an analytic function of 5 also around gc- 
We will come back to this point later. 

Equation ^ gives us the leading order singular part of the free energy. Obviously one can expect 
corrections to such leading behaviour. One kind of corrections, so called analytic corrections, have 
the form: 

esing \t\^'" (1 + ait + flat' + ••■)■ (3) 

Analytic corrections can be simply accounted for by noting that the scaling variable z is 
only approximately z « tL^^". Using a more accurate prescription we can write z = 
t{l + bit + b2t^ + . . .) L^/". Setting this "improved" scaling variable z into Eq. Q we obtain Eq. ^ 
with ai = &i (2 — a) and 02 = (2 — a) (2&2 + 61 (1 — a)) /2 This sort of analytic non-Hnearities 
are a consequence of the fact that the RG equations are only approximately linear at the fixed 
point. 

Besides analytic corrections, we can have corrections to Eq. ^ arising from irrelevant terms. This 



terms contribute to the singular part of the free energy with an additive term of the form [21 1 

B±h{t) 

with A not an integer and h an analytic function of t in zero. Such terms can be accounted for by 
a scaling function of the form 

L-(^+^)-'^/t(t)$i (z). 

In the off-critical regime z ±00, we must have $1 (z) ^ jj_^z''d+c+uj)L' g^^j ^ _ 

Summarizing, the scaling form for the singular part of the free energy, including only the least 
irrelevant operator, is 

esing = L-^'+^^ [$0 (z) + L-'^h (t) $1 (z)] . (4) 
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In the computation of the shift exponents, other terms in the expansion of the energy density may 
play an important role. We write these terms as 

F(g)L-(^+^) + .... (5) 

Where the function F is analytic in g also around gc- Note that the next term in such an expansion 
would be of the form M (g) L~^'^~^''^'^\ with M (g) also analytic. Such a term, close to the critical 
point is equivalent to that arising from an irrelevant operator with w = e. Therefore we do not 
include such terms as they are already included in equation (0]). 

In d — 1, when the theory is conformally invariant (and so C — 1), the origin of a ^-dependent 
term of order L~ like the one in can be justified in the following way. In this case it is 
known that [H, at the critical point, the energy density behaves as e = Coo (ffc) — L~^t:cv/6 + • ■ •. 
Here c is the central charge and v the speed of elementary excitations. Now, the effect of the 
perturbation, beside that of opening a gap, will be that of renormalizing the speed of elementary 
excitations v ^ v{g) giving rise to the function F [g) — —ttcv (g) /6. The presence of such terms 
can also be checked in the XY model via an explicit computation (see section [SH]). The appearance 
of a 5-dependent Casimir term, analogous to our F (g) L~(''+''\ has also been demonstrated for the 
Casimir effect when adding a small mass to the electromagnetic field 0|. 

Armed now with equation Q and possible corrections of the form |(5]) , we can now compute shift 
exponents related to different techniques. The recipe is simple. What we have to do is expand 
the equation determining the finite-size pseudo critical point, in the quasi-critical region L <C ^. It 
suffice to expand the equations up to linear terms in {g — gc). Higher order terms in {g — gc) would 
result in corrections which decay more rapidly with L and so can be safely ignored since we are 
interested in the leading behaviour. 

An important issue in the computation of shift exponents is whether the scaling function <i>o (and/or 
$i) contain a Hnear term or not. If a linear term is present in, say $o, it means that the singular 
part of the energy will not be an even function of t close to the critical point. As a consequence, 
different curves cl (g) at different sizes L will cross close to the critical point. In other words one 
could apply the FSCM equation to the energy itself: Olcl (gl) = 0. However this method gives rise 
to a poor convergence rate, with a shift exponent given by l/i/ as usual. The same shift exponent 
is obtained if one appliesjthe FSCM to the perturbation operator, i.e. solving (g^) = as in 
the original formulation j^. 

Let us consider in detail the equation Olgl (gt) = 0. Up to leading order the solution is given by 

, (d + C)^(F(g,) + $o(0)) ^_i;, 

(i-(d + c).)$;,(o) • 

We note that the leading convergence rate is dictated by the combination {F (gc) + ^'o (0)) being 
non-zero. In the energy density these terms correspond to quantities which scale as L-('*+0 in 
the quasi-critical region. We now present a method which gets rid of such contributions and hence 
allows for a faster convergence rate. 



2.1. The homogeneity condition method 



Since the convergence rate in the FSCM is dictated by {F {gc) + $o (0)) being non-zero, one can 
seek for a method such that the terms of order L~(''+^) cancels out exactly from the expression of 
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cl- This is achieved by requiring that the function Olsl be homogeneous of degree {d + ( + 1). 
We arrive in this way at the homogeneity condition method (HCM), which consists of finding the 
solution to the following equation 



id + C+l)dLeL {9) + LdleL (g) = 0. 



(6) 



Taking into account the full Ansatz Eqs. ([4]) and |(5]) we find that the convergence rate of the 
pseudo-critical sequence so defined, at leading order is given by 



*.HCM 

9l 



9c 



v^u:{d + C + ^)<^i (0)/i(0) 

[(rf + C)^-i]'i>[,(0) 



(7) 



We see that the shift exponent in this case is 9hcm — uo + It should be pointed out that, 
by considering higher derivatives of cl with respect to L, the procedure could be iterated and one 
can in principle find even faster converging sequences. However then one would require explicit 
knowledge of other exponents which are actually unknown. Note also that, despite the origin of 
the HCM is that of eliminating a term of order L~^''+^' from the expression of ei,, the HCM works 
irrespective of the presence of such a term in e^. 

The HCM equation ^ has been written using formally continuous value oi L. To obtain 
the HCM condition at finite size we first replace derivatives with centered finite differences 
which have fast convergence properties, so Dlcl = {cl+sl — cl-sl) I (2(5i) and D\eL = 
(eL+i5L — 2eL + cl-sl) / {5L)^ . Then in principle any combination which reduces to ([6]) for L — > oo 
would produce a pseudo-critical sequence with shift exponent given by Ohcm- However, at finite 
sizes, we can require that the term of order O (L^'^''+''^) in ([5]) cancels out exactly. This choice 
will result in sequences with nicer converging properties. This is achieved by choosing a proper 
combination of DlBl and D\eL , and the HCM equation at finite size becomes 



Pi (i) Dlcl {g) + P2 (L) DlcL (g) - 
where the two polynomials P1/2 are given by 

Pi (i) = \2L'^+'^ (L + SLf+^ + 2L'^+^ (L - 5Lf+'^ - 4 (i^ _ 5L 

P2 (i) = 5L (i + 5Lf+^ - 1"^+^ {L - 5Lf+'^ 

For instance, when d + ( — 2 the HCM equation reduces to — {5L)^ Dl^l {9) + L^Dj^eL (5) 
0, which correctly reduces to |[6)) in the limit of large L. 



(8) 



2.1.1. No linear terms in the scaling functions. As can be inferred by looking at equation ([7]), the 
HCM as is cannot be applied when the scaling function $0 does not contain a linear term. In this 
case one would have $0 (0) = making equation ^ meaningless. However in this situations, finite 
size corrections are smaller, and so it should be possible to envisage a method to turn this point 
into advantage. This is easily achieved by applying the HCM equation ^ to the perturbation 6l 
rather than to the energy density. In this case the solution of (d + C + 1) ^l^l (5) + Ld^DL (5) = 
gives 

*,HCM _ _ (l~M(l-^(rf + C + ^))^(0)1>;(0) , -1/.-. 
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Method 


$0 linear 


*o(0) - 


= 0,$'i(0) = 


FSCM(eL) 


1/^ 


1/v — Lo or not applicable 


not applicable 


FSCM(6l) 


1/^ 


min {2/v, l/v + lu) 


2/v 


HCM(eL) 


l/v + bJ 


\lv 


not applicable 


HCM(6l) 


Ijy 


l/v + uj 


2/v + uj 


MHC 


l/u 


1/v + uj 


2/v + u 



Table 1. Summary of values of shift exponents in different cases. Only the minimum possible 
shift exponent is given. For instance, when also (0) = 0, the shift exponent for HCM(f)i^) 
(which means HCM applied to b^) and MHC becomes at least 2/v + ui. When the method 
formally returns a sequence that does not converge to the critical point in the limit L ^ oo we 
have reported "not applicable". 

and the shift exponent is l/i^ + If also $i does not contain a linear term the shift exponents 
becomes 2/i/ + w when no other cancellations are present. The numerical test we have performed 
(see section [321 below) fall in this case (see also table [J). For completeness we note that, when 
$0 (0) = 0, the FSCM applied to the energy gives formally a shift exponent l/v — to. This is clearly 
a poor convergence and the method can be successfully applied only when \/v > u. 

2.2. The Berezinskii-Kosterlitz-Thouless transition 

The presence of the additional uj in 9hcm suggests that the method could work even in the 
situation where formally v — oo like in the Berezinskii-Kosterlitz-Thouless (BKT) transition. In the 
BKT case the correlation length diverges with an essential singularity of the form ^ ^ exp (at"'"), 
where a and a are constants. For instance in the classical two dimensional XY model a = 1/2. 
Using a modified Ansatz adapted to the BKT, one can show [3| that the HCM Eq. ([H still 
defines a sequence of pseudo-critical points converging to the true critical point with shift exponent 
^Hcli — ^/ {n — 1 + a) with e > and n is an integer greater than one. 

For the application of the HCM equation one simply requires the knowledge of 6^ which is the 
average of a (typically local) quantity over the ground state, and the dynamic exponent C- We 
point out here that, contrary to PRC method (the current method of choice for estimating quantum 
critical points) no need of excited states is required. This is a big advantage since, in every numerical 
method excited states (if available) are typically obtained with smaller accuracy than the ground 
state. Instead the local average can be obtained with very high precision. In addition, the HCM 
is not restricted to quantum mechanical systems and is superior to the PRG in that it produces a 
faster converging sequence. 

2.3. The modified homogeneity conditions 

The dynamic exponent C is one in most known cases, however when not know a priory can be easily 
found by looking at the gap behaviour A — const, x L~'» via exact diagonalization on small lattices. 
However in some cases (like two-dimensional quantum mechanics) assessing C can be a problem, 
but this requirement can in fact be lifted. Considering also the finite size sequence of the ground 
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state energy cl, it is possible to find a modified homogeneity condition (MHC) given by 

bL{gl) - hL^sdgl) _ bL+5Ligl) - ^dgl) 



eUgl) - eL-sLigl) eL+sdgl) - edgl) ' 

or in the continuum version 

dhedg) dleUg)' 



(9) 



(10) 



Using the appropriate scaling Ansatz for the energy cl and the perturbation 6l one can show, that 
the MCH Eq. ^ produces a pseudo-critical sequence g*£^'^'" = gc + const, x L^^mhc _ jf j^^g 
a linear term, the shift exponent is only l/v, growing to 1/u + uj when $q (0) = 0, which becomes 
+ w if also $i (0) = 0. 



2.4- General prescription 

In table [T] we summarize the predicted shift exponents associated to different methods. Notice 
that only the minimal shift exponents are computed. When one has additional cancellations in 
the scaling functions, shift exponents tend to be higher. For example, when neither $o or $i 
contain a Hnear term, the shift exponent computed using the HCM and the MHC is 2/^ + 1^. From 
the available information about the coefficients ^'q{0) and ^"'1(0)1 possibly gathered by means of 
analytical arguments (symmetries, dualities, etc.) and/or graphical inspection of the numerical 
data, one should be able to use the best possible "recipe" to locate the critical point with the 
largest shift exponent. For instance, if $0 carries a linear term in the scaling variable, then the best 
strategy turns out to be the HCM applied to the energy density. 

However, sticking to the HCM, one could also adopt a more "automatic" procedure, in the sense 
that both cl and 6l are computed at once without any previous knowledge about the linear terms 
in the scaling functions. Then, one of the two following situations will occur: (i) Both HCM(ei,) 
and HCM(5i) converge to the same critical point (within numerical uncertainties), meaning that 
^'i (0) ^ and one has to refer to the second or third column of table [H A best fit on the sequence 
with faster convergence then yields the exponent + while the other sequence gives l/v and so 
the two exponents f and uj can be estimated independently, (ii) If instead HCM(eL) and HCM(6l) 
converge to different points (or the energy one does not converge at all), then one has to rely on 
the latter and the best fit returns 2/v + (fifth column of table [1]). Even if now v and uj are no 
longer separated, the convergence towards the critical point is enhanced by a larger exponent. 

This summary concludes our theoretical analysis of the improved shift exponents based on 
the contributions beyond the scaling hypotesis. In the following we will consider practical 
implementations of the accelerated methods HCM and MHC and compare them against previous 
methods namely the PRC and the FSCM. We will consider an exactly solvable case where all the 
quantities can be calculated analytically, a non-integrable example and finally one instance of the 
extremely severe BKT case. In all cases the HCM and MHC represent an improvement with respect 
to previous methods and allow one to reduce to shorter sizes. 
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3. Examples 

3.1. An exactly solvable case: the XY model in transverse field 

To check our methods analytically we consider the one-dimensional XY model in transverse field 
[l^ . It is given by the following Hamiltonian 



H 



L r 

E 

i=i 



V V 



ha' 



(11) 



As usual we will employ PBC to avoid the appearance of surface corrections, and we will use L even. 
The model Eq. ifTTj) can be exactly solved by first passing to spinless fermion via Jordan- Wigner 
transformation and then diagonalizing the quadratic Hamiltonian with a suitable Bogoliubov 
transformation. It turns out that, for all r] ^ 0, there is a critical point at h = he = I which 
belongs to the same universality class as the two-dimensional classical Ising model. All quantities 
of interest like finite size energy density or perturbation 6l, and even lowest gap can be 
computed exactly. In the quasi-critical regime ^ ^ L one is able to find asymptotic expressions 
which are Laurent expansion in L. More details about the calculations can be found in here we 
only quote the results. The energy density at finite size in the quasi-critical region has the following 
form 



6 360 4|r;| 

_2 _ 77r3(3-2 7^ 
12|?/r 2880 177 

(/i-l)Mn(L)+7c + ln(8|77|/7r)-l , ^ / ,^ ^,3 



(h-l) 



+ 0[{h~lf). (12) 



2 7r|7y| 

The Ising model has a = and a logarithmic singularity which at this level shows up in the InL 
term at second order. The logarithmic singularity in principle requires a modified scaling Ansatz 
than that outlined in section [2] (see e.g. [j, [l^). However, stopping at first order, we can still 
compare this expression with our scaling Ansatze Eqns. ^ and ^ expanded in the quasi-critical 
region. We observe that, in this case $q (0) = ^[ (0) = 0. We can also read off the function 
F (h) (here h plays the role of g). Clearly F (h) = -tt |7y| /6 - tt (/i - 1) / (12 \r]\) + O (h ~ if. 
This is in perfect agreement with the conformal field theory (CFT) formula F (h) = —ttcv (h) /6 
taking into account that the central charge is c — 1/2 and the one particle dispersion is 

Afe = 2^/772 + h^ + {l- 772) cos {kf + 2h cos (k) giving rise to an /i-dependent velocity v (h) — 

2 \ri\ + {h — 1) / \r]\ + O {h — if close to the critical point. 

For what regards the terms of order since d + = 2 one deduces that u = 2 (this exponent 
was also noted in |2Q||). Since v — 1 one would conclude that this term arises from an irrelevant 
operator of dimension (2 -I- 2) = 4. A deeper understanding of the finite size corrections of the XY 
model at the critical point {h = 1) can be obtained from the point of view of perturbed CFT. It was 
shown in that the leading finite-size corrections of the critical XY model with PBC have scaling 
dimension 3 (this is an operator belonging to the tower of the energy density in CFT language) and 
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4 (tower of the identity). The operator with scaling dimension 3 would contribute to the energy 
density at finite size with a term of order L~^. However at first order in perturbation theory the 
relevant matrix element vanishes identically, the first non-zero contribution is at second order giving 
rise to a correction |[l8|. Summarizing, the leading corrections arising from irrelevant operators 
with scaHng dimension 3 and 4 are both of order L~* in accordance with Eq. (fT2l) for /i = 1. 

After this short digression we can use Eq. (fT2l) to calculate shift exponents associated to different 
techniques. Applying the HCM conditions one obtains 



72O772 



telling us that 9hcm = 4. Since, from the exact solution one has v = 1 and we just recovered 
w = 2, this result is consistent with 'Ohcm = uj-\-2/v valid when $0 (0) = <f>']^ (0) = 0. The modified 
method MHC gives instead 



720772 



which again confirms Omhc = &hcm = 4 for this case. In this model we can also find exactly the 
shift exponent related to the PRG method. All we need is the expression for the lowest gap at finite 
size, close to the critical point which we recall here from jjJl] 

A,= (.-l)H-^fi2l±A^L- (13) 

H- ■^'P('-l)-2(2 + l.),'+8,']^_, ^ ^ 

192|?7|^ ^ ^ 

Since d = ( = 1 the PRG method consists of imposing Ol [LAl (h*^)] ~ 0, which yields [l3| 

^ * Pi?,G ^ ^ \i — 7^-3 ^ 

^ 48 H 

Finally the FSCM gives h\ = \^ L^'^n'^ in accordance with OpscM = '^/v. 

All in all these explicit calculations confirm that the accelerated HCM and MHC methods are 
superior to the PRG or the FSCM methods. 



3.2. Numerical example: c — 1 transition 

We consider here a spin-1 model with anisotropics given by the following Hamiltonian 

L 

H^J2 [^i ^1+1 + ^J^jVi + + D{S^' f] ' (14) 

i=i 

where as usual, PBC are used. The Hamiltonian Eq. lfT4|) has been used to describe magnetic 
properties of different quasi ID compounds and an extensive literature is available (see e.g. Refs. (lol . 
113|,UJ and reference therein). We fix here to the value Jz = 0.5 where a critical point of order 




Figure 1. (color on-line) The sequences of pseudo-critical points obtained with the FSCM, the 
PRG and new methods HCM (given by Eq. Jsj) and the MHC Eq. JgJ. Continuous lines are 
algebraic best-fits to the data. 





FSCM 


PRG 


HCM 


MHC 




0.647 


0.640 


0.630 


0.633 


fitted e 


0.79 


1.50 


7.6 


7.6 



Table 2. Extrapolated critical points together with shift exponents for various methods. 



Dc = 0.6 0.7 has been estimated previously. At this particular transition point it has been 
previously found a pretty large value of the correlation length exponent v = 2.38, which implies 
a slow scaling in l/v. 

The advantage of having fast converging sequences is that one can limit the numerical calculations 
to short lattice sizes that can be reached by exact diagonalization. To illustrate our procedures we 
used exact Lanczos diagonalizations for chain of length L = 8, 10, 12, 14, 16 and used high precision 
DMRG simulations (with as many as 3'' optimized states and many zips) only for chain length 
L = 18, 20. The role of 6l is given here by 6l = {{Sff). We tested different methods: the FSCM, 
PRG, HCM, and MHC. We first observed that the energy curves at different size do not cross. 
This means that there are no linear terms in the scaling functions and consequently we applied the 
FSCM and the HCM to the perturbation operator 6/,. The results are plotted in Figure [TJ It is 
apparent that the HCM and MHC result in much more rapidly converging sequences respect to 
other methods, included the PRG. 

The extrapolated critical points together with the relative shift exponents for different methods, 
are summarized in Table [2l It is apparent that HCM and MHC produce a very fast converging 
sequence with similar shift exponents given by Ohcm ~ Omhc ~ 7.6. 
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Figure 2. (color on-line) Left panel: Search for the zeros of Eq. |[8j with 5L = 4. Right Panel: 
Extrapolations of the pseudo-critical points J| (cyan line for all the points, blue line with L = 12 
excluded); the black cross marks the accepted value 0.2411 (see text for further details). 



3.3. BKT transition 

To finish our series of numerical test we consider the example of the extremely severe case of BKT 
transition. In such a transition the gap opens up so slowly (remind A ~ ~ exp {—at~'^)) that 
the PRG method is practically not applicable. As an example of BKT transition we consider the 
spin-1/2 Heisenberg model with frustration due to next-to-nearest neighbors interaction: 

L 

The model is equivalent to a 2-legs zigzag ladder with L/2 rungs. We will also fix Ji = 1 for 
clarity. The model is gapless for J2 < J2C while it has a doubly degenerate ground state for 
J2 > J2c- The best estimate for the critical point was done by Okamoto and Nomura in [16] using 
a model-specific crossing method which requires the computations of highly excited states. They 
found J2c — 0.2411 ±0.0001. The role of the parameter g is taken here by J2 while the perturbation 
operator is given by — {Sj-Sj+2). The numerical calculations were done using ladders with PBC 
and up to L/2 — 16 rungs, an effective SL — 4, and 1024 DMRG states. We now pass to examine 
the different methods. As already mentioned the PRG is not feasible because of the extremely slow 
opening of the gap. The FSCM also fails as it turns out to be impossible to find a real solution to 
the FSCM equation. Instead the HCM results in a nice pseudo-critical sequence converging with a 
shift exponent close to one. The resulting sequence is plotted in Figured A linear fit including all 
points from L = 12 to L = 32 gives J2C — 0.242 ± 0.006 while excluding the point at L = 12 we 
obtain a more precise value J2C — 0.2553 ± 0.0008 although larger than the one reported in [l^ . 
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4. Conclusions 

Locating the critical point is one of the first issues which arise when considering a statistical 
mechanical problem. Not much attention has been devoted so far in searching for efficient 
methods to extract critical points from numerical simulations. Available methods are, among 
others: i) looking for the maxima of a susceptibility ii) intersection points of reduced (Binder's) 
cumulants at different system sizes, iii) look for closing of gaps in quantum systems (PRG). These 
methods produce sequences of pseudo-critical points which converge to the true critical point as 
9l = Sc + O (L^^) with a leading algebraic term controlled by the shift exponent 6. The general 
behef is that the shift exponent is given by the inverse of the correlation length exponent: 9 = l/v. 
In this article we developed accelerated methods for the location of critical points. Our methods 
can be applied at any second order transition in any dimension quantum or classical, temperature 
or parameter driven. The resulting pseudo-critical sequences are characterized by a large value of 
the shift exponent which enables these methods to be efficient even in the extremely hard case of 
BKT transition where formally 'V = cxd". 

The advantage of having such methods at disposal, is that one can limit the analysis to much 
shorter system sizes. This is particularly evident in two dimensional quantum mechanical problems 
with frustration for which quantum Monte Carlo methods are plagued with the sign problem and 
the reliable numerical data are therefore restricted to very short sizes. 
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